home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 1.iso / toolbox / src / demos / OpenGL / backtrace / Unitdisk.c++ < prev    next >
C/C++ Source or Header  |  1996-11-11  |  21KB  |  867 lines
  1. /*
  2.  * (c) Copyright 1993, Silicon Graphics, Inc.
  3.  * ALL RIGHTS RESERVED 
  4.  * Permission to use, copy, modify, and distribute this software for 
  5.  * any purpose and without fee is hereby granted, provided that the above
  6.  * copyright notice appear in all copies and that both the copyright notice
  7.  * and this permission notice appear in supporting documentation, and that 
  8.  * the name of Silicon Graphics, Inc. not be used in advertising
  9.  * or publicity pertaining to distribution of the software without specific,
  10.  * written prior permission. 
  11.  *
  12.  * THE MATERIAL EMBODIED ON THIS SOFTWARE IS PROVIDED TO YOU "AS-IS"
  13.  * AND WITHOUT WARRANTY OF ANY KIND, EXPRESS, IMPLIED OR OTHERWISE,
  14.  * INCLUDING WITHOUT LIMITATION, ANY WARRANTY OF MERCHANTABILITY OR
  15.  * FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT SHALL SILICON
  16.  * GRAPHICS, INC.  BE LIABLE TO YOU OR ANYONE ELSE FOR ANY DIRECT,
  17.  * SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY
  18.  * KIND, OR ANY DAMAGES WHATSOEVER, INCLUDING WITHOUT LIMITATION,
  19.  * LOSS OF PROFIT, LOSS OF USE, SAVINGS OR REVENUE, OR THE CLAIMS OF
  20.  * THIRD PARTIES, WHETHER OR NOT SILICON GRAPHICS, INC.  HAS BEEN
  21.  * ADVISED OF THE POSSIBILITY OF SUCH LOSS, HOWEVER CAUSED AND ON
  22.  * ANY THEORY OF LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE
  23.  * POSSESSION, USE OR PERFORMANCE OF THIS SOFTWARE.
  24.  * 
  25.  * US Government Users Restricted Rights 
  26.  * Use, duplication, or disclosure by the Government is subject to
  27.  * restrictions set forth in FAR 52.227.19(c)(2) or subparagraph
  28.  * (c)(1)(ii) of the Rights in Technical Data and Computer Software
  29.  * clause at DFARS 252.227-7013 and/or in similar or successor
  30.  * clauses in the FAR or the DOD or NASA FAR Supplement.
  31.  * Unpublished-- rights reserved under the copyright laws of the
  32.  * United States.  Contractor/manufacturer is Silicon Graphics,
  33.  * Inc., 2011 N.  Shoreline Blvd., Mountain View, CA 94039-7311.
  34.  *
  35.  * OpenGL(TM) is a trademark of Silicon Graphics, Inc.
  36.  */
  37. #include <GL/glu.h>
  38. #include <stdio.h>
  39. #include <stdlib.h>
  40. #include <math.h>
  41.  
  42. #include "Unitdisk.h"
  43.  
  44. const float fudge = .000001;
  45. const float M_2PI = 2.0 * M_PI;
  46.  
  47. Unitdisk::Unitdisk()
  48. {
  49.   rdivisions = 3;
  50.   tdivisions = 10;
  51.   points = normals = NULL;
  52.   colors = NULL;
  53.   points_size = normals_size = colors_size = 0;
  54.   angle = M_PI;
  55.   zaxis[0] = 0;
  56.   zaxis[1] = 0;
  57.   zaxis[2] = 1;
  58.  
  59.   still_in_xy = 1;
  60.  
  61.   sintable = costable = NULL;
  62. }
  63.  
  64. Unitdisk::~Unitdisk()
  65. {
  66. }
  67.  
  68. void Unitdisk::draw()
  69. {
  70.   if (points == NULL) return;
  71.   glNormal3f(0, 0, 1);
  72.   if (colors == NULL) draw_nocolors();
  73.   else if (normals == NULL) draw_colors_nonormals();
  74.   else draw_colors_normals();
  75. }
  76.  
  77. void Unitdisk::draw_nocolors()
  78. {
  79.   int r, t, p1, p2;
  80.   int has_n;
  81.  
  82.   has_n = (normals != NULL);
  83.  
  84.   for (t = 1; t < tdivisions; t++) {
  85.     glBegin(GL_QUAD_STRIP);
  86.     p1 = t * (rdivisions + 1);
  87.     p2 = (t - 1) * (rdivisions + 1);
  88.     for (r = 0; r <= rdivisions; r++, p1++, p2++) {
  89.       if (has_n) glNormal3fv(normals[p1].pt);
  90. #ifdef TEXTURE
  91.       glTexCoord2fv(points[p1].pt);
  92. #endif
  93.       glVertex3fv(points[p1].pt);
  94.       if (has_n) glNormal3fv(normals[p2].pt);
  95. #ifdef TEXTURE
  96.       glTexCoord2fv(points[p2].pt);
  97. #endif
  98.       glVertex3fv(points[p2].pt);
  99.     } 
  100.     glEnd(); 
  101.   }
  102.   
  103.   glBegin(GL_QUAD_STRIP); 
  104.   p1 = 0;
  105.   p2 = (rdivisions + 1) * (tdivisions - 1);
  106.   for (r = 0; r <= rdivisions; r++, p1++, p2++) {
  107.     if (has_n) glNormal3fv(normals[p1].pt);
  108. #ifdef TEXTURE
  109.     glTexCoord2fv(points[p1].pt);
  110. #endif
  111.     glVertex3fv(points[p1].pt);
  112.     if (has_n) glNormal3fv(normals[p2].pt);
  113. #ifdef TEXTURE
  114.     glTexCoord2fv(points[p2].pt);
  115. #endif
  116.     glVertex3fv(points[p2].pt);
  117.   }
  118.   glEnd();
  119. }
  120.  
  121. void Unitdisk::draw_colors_nonormals()
  122. {  
  123.   int r, t, p1, p2;
  124.  
  125.   for (t = 1; t < tdivisions; t++) {
  126.     glBegin(GL_QUAD_STRIP);
  127.     p1 = t * (rdivisions + 1);
  128.     p2 = (t - 1) * (rdivisions + 1);
  129.     for (r = 0; r <= rdivisions; r++, p1++, p2++) {
  130.       glColor4fv(colors[p1].c);
  131. #ifdef TEXTURE
  132.       glTexCoord2fv(points[p1].pt);
  133. #endif
  134.       glVertex3fv(points[p1].pt);
  135.       glColor4fv(colors[p2].c);
  136. #ifdef TEXTURE
  137.       glTexCoord2fv(points[p2].pt);
  138. #endif
  139.       glVertex3fv(points[p2].pt);
  140.     } 
  141.     glEnd(); 
  142.   }
  143.   
  144.   glBegin(GL_QUAD_STRIP); 
  145.   p1 = 0;
  146.   p2 = (rdivisions + 1) * (tdivisions - 1);
  147.   for (r = 0; r <= rdivisions; r++, p1++, p2++) {
  148.     glColor4fv(colors[p1].c);
  149. #ifdef TEXTURE
  150.     glTexCoord2fv(points[p1].pt);
  151. #endif
  152.     glVertex3fv(points[p1].pt);
  153.     glColor4fv(colors[p2].c);
  154. #ifdef TEXTURE
  155.     glTexCoord2fv(points[p2].pt);
  156. #endif
  157.     glVertex3fv(points[p2].pt);
  158.   }
  159.   glEnd();
  160. }
  161.  
  162. void Unitdisk::draw_colors_normals()
  163. {
  164.   int r, t, p1, p2;
  165.  
  166.   for (t = 1; t < tdivisions; t++) {
  167.     glBegin(GL_QUAD_STRIP);
  168.     p1 = t * (rdivisions + 1);
  169.     p2 = (t - 1) * (rdivisions + 1);
  170.     for (r = 0; r <= rdivisions; r++, p1++, p2++) {
  171.       glColor4fv(colors[p1].c);
  172.       glNormal3fv(normals[p1].pt);
  173. #ifdef TEXTURE
  174.       glTexCoord2fv(points[p1].pt);
  175. #endif
  176.       glVertex3fv(points[p1].pt);
  177.       glColor4fv(colors[p2].c);
  178.       glNormal3fv(normals[p2].pt);
  179. #ifdef TEXTURE
  180.       glTexCoord2fv(points[p2].pt);
  181. #endif
  182.       glVertex3fv(points[p2].pt);
  183.     } 
  184.     glEnd(); 
  185.   }
  186.   
  187.   glBegin(GL_QUAD_STRIP); 
  188.   p1 = 0;
  189.   p2 = (rdivisions + 1) * (tdivisions - 1);
  190.   for (r = 0; r <= rdivisions; r++, p1++, p2++) {
  191.     glColor4fv(colors[p1].c);
  192.     glNormal3fv(normals[p1].pt);
  193. #ifdef TEXTURE
  194.     glTexCoord2fv(points[p1].pt);
  195. #endif
  196.     glVertex3fv(points[p1].pt);
  197.     glColor4fv(colors[p2].c);
  198.     glNormal3fv(normals[p2].pt);
  199. #ifdef TEXTURE
  200.     glTexCoord2fv(points[p2].pt);
  201. #endif
  202.     glVertex3fv(points[p2].pt);
  203.   }
  204.   glEnd();
  205. }
  206.  
  207. void Unitdisk::draw_by_perimeter(int pass_colors, int pass_norms, 
  208.                  int pass_tex)
  209. {
  210.   int i, index, r = rdivisions + 1;
  211.  
  212.   if (points == NULL) return;
  213.   if (pass_colors && colors == NULL) {
  214.     fprintf(stderr, "Warning:  No colors to draw in Unitdisk.c++");
  215.     pass_colors = 0;
  216.   }
  217.   if (pass_norms && normals == NULL) {
  218.     fprintf(stderr, "Warning:  No normals to draw in Unitdisk.c++");
  219.     pass_norms = 0;
  220.   }
  221.   glBegin(GL_POLYGON);
  222.   for (i = 0, index = rdivisions; i < tdivisions; i++, index += r) {
  223.     if (pass_colors) glColor4fv(colors[index].c);
  224.     if (pass_norms) glNormal3fv(normals[index].pt);
  225. #ifdef TEXTURE
  226.     if (pass_tex) glTexCoord2fv(points[index].pt);
  227. #endif
  228.     glVertex3fv(points[index].pt);
  229.   }
  230.   glEnd();
  231.   
  232. }
  233.  
  234. void Unitdisk::set_angle(float new_angle)
  235. {
  236.   angle = new_angle;
  237. }
  238.  
  239. GLfloat Unitdisk::get_angle()
  240. {
  241.   return angle;
  242. }
  243.  
  244. GLfloat Unitdisk::get_radius()
  245. {
  246.   return cos((M_PI - angle) / 2.0);
  247. }
  248.  
  249. void Unitdisk::set_divisions(int new_rdivisions, int new_tdivisions)
  250. {
  251.   if (tdivisions != new_tdivisions) {
  252.     delete sintable;
  253.     delete costable;
  254.     sintable = costable = NULL;
  255.   }
  256.   if (tdivisions != new_tdivisions || rdivisions != new_rdivisions) {
  257.     rdivisions = new_rdivisions;
  258.     tdivisions = new_tdivisions;
  259.     free_points();
  260.     free_normals();
  261.     free_colors();
  262.   }
  263. }
  264.   
  265. int Unitdisk::get_rdivisions()
  266. {
  267.   return rdivisions;
  268. }
  269.  
  270. int Unitdisk::get_tdivisions()
  271. {
  272.   return tdivisions;
  273. }
  274.  
  275. void Unitdisk::alloc_points()
  276. {
  277.   int npoints = get_npoints();
  278.   if (npoints > points_size) {
  279.     delete points;
  280.     points = new Point[npoints];
  281.     points_size = npoints;
  282.   }
  283. }
  284.  
  285. void Unitdisk::alloc_normals()
  286. {
  287.   int npoints = get_npoints();
  288.   if (npoints > normals_size) {
  289.     delete normals;
  290.     normals = new Point[npoints];
  291.     normals_size = npoints;
  292.   }
  293. }
  294.  
  295. void Unitdisk::alloc_points_normals()
  296. {
  297.   alloc_points();
  298.   alloc_normals();
  299. }
  300.  
  301. void Unitdisk::free_points() 
  302. {
  303.   delete points;
  304.   points = NULL;
  305.   points_size = 0;
  306. }
  307.  
  308. void Unitdisk::free_normals()
  309. {
  310.   delete normals;
  311.   normals = NULL;
  312.   normals_size = 0;
  313. }
  314.  
  315. void Unitdisk::free_points_normals()
  316. {
  317.   free_points();
  318.   free_normals();
  319. }
  320.  
  321. void Unitdisk::fill_points()
  322. {
  323.   alloc_points();
  324.   fill_either(points);
  325. }
  326.  
  327. void Unitdisk::fill_normals()
  328. {
  329.   alloc_normals();
  330.   fill_either(normals);
  331. }
  332.  
  333. void Unitdisk::fill_points_normals()
  334. {
  335.   alloc_points();
  336.   alloc_normals();
  337.   fill_either(points);
  338.   fill_either(normals);
  339. }
  340.  
  341. void Unitdisk::copy_points(Unitdisk src) 
  342. {
  343.   set_divisions(src.rdivisions, src.tdivisions);
  344.   alloc_points();
  345.   copy_either(points, src.points);
  346. }
  347.  
  348. void Unitdisk::copy_normals(Unitdisk src) 
  349. {
  350.   set_divisions(src.rdivisions, src.tdivisions);
  351.   alloc_normals();
  352.   copy_either(normals, src.normals);
  353. }
  354.  
  355. void Unitdisk::copy_either(Point *dpt, Point *spt) {
  356.   int i, npoints;
  357.   npoints = get_npoints();
  358.   for (i = 0; i < npoints; i++, dpt++, spt++) {
  359.     dpt->pt[0] = spt->pt[0];
  360.     dpt->pt[1] = spt->pt[1];
  361.     dpt->pt[2] = spt->pt[2];
  362.   }
  363. }  
  364.  
  365. void Unitdisk::copy_normals_from_points(Unitdisk src) {
  366.   set_divisions(src.rdivisions, src.tdivisions);
  367.   alloc_normals();
  368.   copy_either(normals, src.points);
  369. }
  370.  
  371. void Unitdisk::copy_normals_from_points() {
  372.   copy_normals_from_points(*this);
  373. }
  374.  
  375. void Unitdisk::fill_either(Point *what) {
  376.   int t, r;
  377.   int i;
  378.  
  379.   fill_either_strip1(what);
  380.   if (sintable == NULL) fill_trig_tables();
  381.   i = rdivisions + 1;
  382.   for (t = 1; t < tdivisions; t++) {
  383.     for (r = 0; r <= rdivisions; r++) {
  384.       what[i].pt[0] = costable[t] * what[r].pt[0];
  385.       what[i].pt[1] = sintable[t] * what[r].pt[0];
  386.       what[i].pt[2] = what[r].pt[2];
  387.       i++;
  388.     }
  389.   }
  390. }
  391.  
  392. void Unitdisk::fill_points_strip1()
  393. {
  394.   alloc_points();
  395.   fill_either_strip1(points);
  396. }
  397.  
  398. void Unitdisk::fill_either_strip1(Point *what) {
  399.   float radius, rinc;
  400.   int r;
  401.  
  402.   rinc = get_radius() / (float)rdivisions;
  403.   radius = 0.0;
  404.   for (r = 0; r <= rdivisions; r++, radius += rinc) {
  405.     what[r].pt[0] = radius;
  406.     what[r].pt[1] = 0;
  407.     // Round-off error avoidance hack
  408.     what[r].pt[2] = 1.0 - what[r].pt[0]*what[r].pt[0];
  409.     if (what[r].pt[2] > 0.0) what[r].pt[2] = sqrt(what[r].pt[2]);
  410.     else what[r].pt[2] = 0.0;
  411.  }   
  412. }
  413.  
  414. void Unitdisk::translate(Point trans)
  415. {
  416.   int i, npoints;
  417.   npoints = get_npoints();
  418.   for (i = 0; i < npoints; i++) {
  419.     points[i].pt[0] += trans.pt[0];
  420.     points[i].pt[1] += trans.pt[1];
  421.     points[i].pt[2] += trans.pt[2];
  422.   }
  423. }
  424.  
  425. void Unitdisk::scale(float s) 
  426. {
  427.   int i, npoints;
  428.   npoints = get_npoints();
  429.   for (i = 0; i < npoints; i++) {
  430.     points[i].pt[0] *= s;
  431.     points[i].pt[1] *= s;
  432.     points[i].pt[2] *= s;
  433.   }
  434. }
  435.  
  436. void Unitdisk::scale_translate(float s, Point trans) 
  437. {
  438.   int i, npoints;
  439.   npoints = get_npoints();
  440.   for (i = 0; i < npoints; i++) {
  441.     points[i].pt[0] = points[i].pt[0]*s + trans.pt[0];
  442.     points[i].pt[1] = points[i].pt[1]*s + trans.pt[1];
  443.     points[i].pt[2] = points[i].pt[2]*s + trans.pt[2];
  444.   }
  445. }
  446.  
  447. int Unitdisk::get_npoints() 
  448. {
  449.   return (rdivisions + 1) * tdivisions;
  450. }
  451.  
  452. void Unitdisk::project()
  453. {
  454.   int i, npoints;
  455.  
  456.   if (normals == NULL) {
  457.     fprintf(stderr, "Warning:  No normals defined when project() called.\n");
  458.     fill_normals();
  459.   }
  460.   if (points == NULL) fill_points();
  461.   npoints = get_npoints();
  462.   for (i = 0; i < npoints; i++) {
  463.     /* I'm not sure quite what the justification for this is, but it
  464.      * seems to work */
  465.     if (normals[i].pt[2] < 0.0) normals[i].pt[2] = -normals[i].pt[2];
  466.     points[i] = points[i].project_direction(normals[i]);
  467.   }
  468. }
  469.  
  470. void Unitdisk::project(Point projpt)
  471. {
  472.   int i, npoints = get_npoints();
  473.   float x, y, z;
  474.   Point *pt;
  475.  
  476.   if (points == NULL) fill_points();
  477.   x = projpt.pt[0];
  478.   y = projpt.pt[1];
  479.   z = projpt.pt[2];
  480.   pt = points;
  481.   for (i = 0; i < npoints; i++, pt++) pt->project_self(x, y, z);
  482. }
  483.  
  484. void Unitdisk::project_borrow_points(Unitdisk src) 
  485. {
  486.   int i, npoints = get_npoints();
  487.   Point *pt, *spt, *sn;
  488.   spt = src.points;
  489.   sn = normals;
  490.   alloc_points();
  491.   pt = points;
  492.   for (i = 0; i < npoints; i++, pt++, spt++, sn++) {
  493.     /* I'm not sure quite what the justification for this is, but it
  494.      * seems to work */
  495.     if (normals[i].pt[2] < 0.0) normals[i].pt[2] = -normals[i].pt[2];
  496.     pt->compute_projected(spt->pt[0], spt->pt[1], spt->pt[2], 
  497.               sn->pt[0], sn->pt[1], sn->pt[2]);
  498.   }
  499. }
  500.  
  501. void Unitdisk::refract_normals(Point light, GLfloat I)
  502. {
  503.   Point dlight;
  504.   float cos1, sin1, cos2, sin2;
  505.   int use_normals;
  506.   int r, t, i;
  507.  
  508.   if (points == NULL) {
  509.     fprintf(stderr, "Attempting to refract without points.\n");
  510.     fill_normals();
  511.   }
  512.  
  513.   use_normals = (normals != NULL);
  514.   alloc_normals();
  515.  
  516.   /* Do the theta = 0 diagonal */
  517.   for (r = 0; r <= rdivisions; r++) {
  518.     /* Find the original normal - use the unit of the points if there are
  519.      * no normals */
  520.     if (!use_normals) normals[r] = points[r].unit();
  521.  
  522.     /* Compute the direction to the light */
  523.     dlight = (light - points[r]).unit();
  524.  
  525.     /* Compute the cosine and the sine of the original angle */
  526.     cos1 = dlight.dot(normals[r]);
  527.     sin1 = 1.0 - cos1*cos1;
  528.     if (sin1 <= 0.0) continue;
  529.     sin1 = sqrt(sin1);
  530.  
  531.     /* Compute the cosine and the sine of the new angle */
  532.     sin2 = sin1 / I;
  533.     cos2 = sqrt(1.0 - sin2*sin2);
  534.  
  535.     /* Rotate the normal by the new sine and cosine */
  536.     normals[r] = normals[r].rotate_abouty(cos2, -sin2);
  537.   }
  538.  
  539.   /* Copy the rest of the rows from the current row */
  540.   i = rdivisions + 1;
  541.   if (sintable == NULL) fill_trig_tables();
  542.   for (t = 1; t < tdivisions; t++) {
  543.     for (r = 0; r <= rdivisions; r++) {
  544.       normals[i].pt[0] = costable[t] * normals[r].pt[0];
  545.       normals[i].pt[1] = sintable[t] * normals[r].pt[0];
  546.       normals[i].pt[2] = normals[r].pt[2];
  547.       i++;
  548.     }
  549.   }
  550.  
  551. }
  552.  
  553. void Unitdisk::face_direction(Point d)
  554. {
  555.   face_direction(d, *this);
  556. }
  557.  
  558. void Unitdisk::face_direction(Point d, Unitdisk src) 
  559. {
  560.   Point *spt, *dpt, *sn, *dn;
  561.   float sin1, cos1;
  562.   float x;
  563.   int npoints, i;
  564.  
  565.   if (d.pt[1]) {
  566.     fprintf(stderr, 
  567.         "Internal error:  Can't face in direction not in xz plane.");
  568.     return;
  569.   }
  570.  
  571.   cos1 = d.pt[2];
  572.   sin1 = d.pt[0];
  573.   
  574.   if (sin1 * sin1 > fudge) {
  575.     spt = src.points;
  576.     dpt = points;
  577.     sn = src.normals;
  578.     dn = normals;
  579.     /* Change this to be seperate loops for points&&normals, points, normals
  580.      * (faster than testing every iteration */
  581.     npoints = get_npoints();
  582.     for (i = 0; i < npoints; i++) {
  583.       if (points != NULL) {
  584.     x = spt->pt[0];
  585.     dpt->pt[0] = x*cos1 + spt->pt[2]*sin1;
  586.     dpt->pt[1] = spt->pt[1];
  587.     dpt->pt[2] = -x*sin1 + spt->pt[2]*cos1;
  588.     spt++; dpt++;
  589.       }
  590.       if (normals != NULL) {
  591.     x = sn->pt[0];
  592.     dn->pt[0] = x*cos1 + sn->pt[2]*sin1;
  593.     dn->pt[1] = sn->pt[1];
  594.     dn->pt[2] = -x*sin1 + sn->pt[2]*cos1;
  595.     sn++; dn++;
  596.       }
  597.     }
  598.   } else if (points != NULL && points != src.points) {
  599.     copy_points(src);
  600.     if (normals != NULL) copy_normals(src);
  601.   }
  602. }  
  603.  
  604. void Unitdisk::alloc_colors()
  605. {
  606.   int ncolors = get_npoints();
  607.   if (ncolors > colors_size) {
  608.     delete colors;
  609.     colors = new Color[ncolors];
  610.     colors_size = ncolors;
  611.   }
  612. }
  613.  
  614. void Unitdisk::map_normals_to_colors() 
  615. {
  616.   int t, r;
  617.   int i;
  618.  
  619.   if (normals == NULL) fill_normals();
  620.   alloc_colors();
  621.   i = 0;
  622.   for (t = 1; t <= tdivisions; t++) {
  623.     for (r = 0; r <= rdivisions; r++) {
  624.       colors[i] *= normals[r].pt[2];
  625.       colors[i].c[3] = 1;
  626.       i++;
  627.     }
  628.   }
  629. }
  630.  
  631. void Unitdisk::map_z_to_colors()
  632. {
  633.   int i, npoints = get_npoints();
  634.  
  635.   if (points == NULL) {
  636.     fprintf(stderr, "Warning:  no points defined in map_z_to_colors()\n");
  637.     fill_points();
  638.   }
  639.  
  640.   alloc_colors();
  641.   for (i = 0; i < npoints; i++) {
  642.     colors[i] = points[i].pt[2];
  643.     colors[i].c[3] = 1;
  644.   }
  645. }
  646.  
  647. void Unitdisk::scale_alpha_by_z()
  648. {
  649.   int i, npoints = get_npoints();
  650.   
  651.   if (colors == NULL) alloc_colors();
  652.   if (points == NULL) {
  653.     alloc_points();
  654.     fill_points();
  655.   }
  656.   for (i = 0; i < npoints; i++) colors[i].c[3] *= points[i].pt[2];
  657. }
  658.  
  659. void Unitdisk::scale_colors_by_z() 
  660. {
  661.   int i, npoints = get_npoints();
  662.  
  663.   if (colors == NULL) alloc_colors();
  664.   if (points == NULL) {
  665.     alloc_points();
  666.     fill_points();
  667.   }
  668.   for (i = 0; i < npoints; i++) colors[i] *= points[i].pt[2];
  669. }
  670.  
  671. void Unitdisk::scale_colors_by_normals(Point light)
  672. {
  673.   scale_colors_by_normals(light, *this);
  674. }
  675.  
  676. void Unitdisk::scale_colors_by_normals(Point light, Unitdisk src_normals)
  677. {
  678.   scale_colors_by_either(light, src_normals.normals);
  679. }
  680.  
  681. void Unitdisk::scale_colors_by_points(Point light, Unitdisk src_points)
  682. {
  683.   scale_colors_by_either(light, src_points.points);
  684. }
  685.  
  686. void Unitdisk::scale_colors_by_either(Point light, Point *what)
  687. {
  688.   int t, r;
  689.   int i;
  690.   if (what == NULL) {
  691.     fprintf(stderr, "Scaling colors to NULL pointer.\n");
  692.     return;
  693.   }
  694.   if (light.pt[0] || light.pt[1] || light.pt[2] < 0.0) {
  695.     fprintf(stderr, "Light not on z axis in scale_colors_by_normals.\n");
  696.   }
  697.  
  698.   alloc_colors();
  699.   for (r = 0; r <= rdivisions; r++) 
  700.     colors[r] *= what[r].dot(light);
  701.   i = rdivisions + 1;
  702.   for (t = 1; t < tdivisions; t++) {
  703.     for (r = 0; r <= rdivisions; r++) {
  704.       colors[i] = colors[r];
  705.       i++;
  706.     }
  707.   }
  708. }  
  709.  
  710. void Unitdisk::set_colors(Color c)
  711. {
  712.   Color *dst;
  713.   int i, npoints = get_npoints();
  714.   alloc_colors();
  715.   dst = colors;
  716.   for (i = 0; i < npoints; i++, dst++) {
  717.     dst->c[0] = c.c[0];
  718.     dst->c[1] = c.c[1];
  719.     dst->c[2] = c.c[2];
  720.     dst->c[3] = c.c[3];
  721.   }
  722. }
  723.  
  724. void Unitdisk::add_colors(Color c)
  725. {
  726.   int i, npoints = get_npoints();
  727.   alloc_colors();
  728.   for (i = 0; i < npoints; i++) colors[i] += c;
  729. }
  730.  
  731. void Unitdisk::free_colors()
  732. {
  733.   delete colors;
  734.   colors = NULL;
  735.   colors_size = 0;
  736. }
  737.  
  738. inline float Unitdisk::area_triangle(Point a, Point b, Point c)
  739. {
  740.   return (((a.pt[0]*b.pt[1] + b.pt[0]*c.pt[1] + c.pt[0]*a.pt[1]) -
  741.        (a.pt[1]*b.pt[0] + b.pt[1]*c.pt[0] + c.pt[1]*a.pt[0])) * .5);
  742. }
  743.  
  744. inline float Unitdisk::area_triangle(GLfloat *a, GLfloat *b, 
  745.                      GLfloat *c) 
  746. {
  747.   return (((a[0]*b[1] + b[0]*c[1] + c[0]*a[1]) -
  748.        (a[1]*b[0] + b[1]*c[0] + c[1]*a[0])) * .5);
  749. }
  750.  
  751. inline float Unitdisk::area_2triangle(GLfloat *a, GLfloat *b, 
  752.                      GLfloat *c) 
  753. {
  754.   return ((a[0]*b[1] + b[0]*c[1] + c[0]*a[1]) -
  755.       (a[1]*b[0] + b[1]*c[0] + c[1]*a[0]));
  756. }
  757.  
  758. void Unitdisk::scale_colors_by_darea(Unitdisk disk) 
  759. {
  760.   int pt1, pt2, pt3;
  761.   int t, r, i;
  762.   int npoints = get_npoints();
  763.   float *rproducts1, *tproducts1, *rproducts2, *tproducts2;
  764.   float area1, area2;
  765.  
  766.   rproducts1 = new float[npoints];
  767.   tproducts1 = new float[npoints];
  768.   rproducts2 = new float[npoints];
  769.   tproducts2 = new float[npoints];
  770.  
  771.   /* Compute the products of the segments which make up the disk - 
  772.    * these will later be used in the area calculations */
  773.   i = 0;
  774.   for (t = 0; t < tdivisions; t++) {
  775.     for (r = 0; r < rdivisions; r++) {
  776.       pt1 = i;
  777.       pt2 = i + 1;
  778.       rproducts1[i] = (points[pt1].pt[0]*points[pt2].pt[1] -
  779.                points[pt1].pt[1]*points[pt2].pt[0]);
  780.       rproducts2[i] = (disk.points[pt1].pt[0]*disk.points[pt2].pt[1] - 
  781.                disk.points[pt1].pt[1]*disk.points[pt2].pt[0]);
  782.       pt2 = ((t+1)%tdivisions)*(rdivisions + 1) + r;
  783.       tproducts1[i] = (points[pt1].pt[0]*points[pt2].pt[1] -
  784.                points[pt1].pt[1]*points[pt2].pt[0]);
  785.       tproducts2[i] = (disk.points[pt1].pt[0]*disk.points[pt2].pt[1] - 
  786.                disk.points[pt1].pt[1]*disk.points[pt2].pt[0]);
  787.       i++;
  788.     }
  789.     pt1 = i;
  790.     pt2 = ((t+1)%tdivisions)*(rdivisions + 1) + r;
  791.     tproducts1[i] = (points[pt1].pt[0]*points[pt2].pt[1] -
  792.              points[pt1].pt[1]*points[pt2].pt[0]);
  793.     tproducts2[i] = (disk.points[pt1].pt[0]*disk.points[pt2].pt[1] - 
  794.              disk.points[pt1].pt[1]*disk.points[pt2].pt[0]);
  795.     i++;
  796.   }
  797.  
  798.   /* Compute the area at the center of the disk */
  799.   area1 = area2 = 0.0;
  800.   r = 1;
  801.   for (t = 0; t <= tdivisions; t++) {
  802.     pt1 = (t%tdivisions)*(rdivisions+1) + r;
  803.     area1 += tproducts1[pt1];
  804.     area2 += tproducts2[pt1];
  805.   }
  806.   if (area1 != 0.0) area1 = fabs(area2 / area1);
  807.   for (t = 0; t < tdivisions; t++) {
  808.     colors[t*(rdivisions+1)] *= area1;
  809.   }
  810.  
  811.   for (t = 0; t < tdivisions; t++) {
  812.     for (r = 1; r < rdivisions; r++) {
  813.       pt1 = (t ? t-1 : tdivisions-1)*(rdivisions+1) + r - 1;
  814.       pt3 = t*(rdivisions + 1) + r - 1;
  815.       pt2 = ((t+1) % tdivisions)*(rdivisions+1) + r - 1;
  816.       area1 = rproducts1[pt1] + rproducts1[pt1 + 1];
  817.       area1 += tproducts1[pt1 + 2] + tproducts1[pt3 + 2];
  818.       area1 -= rproducts1[pt2 + 1] + rproducts1[pt2];
  819.       area1 -= tproducts1[pt3] + tproducts1[pt1];
  820.       area2 = rproducts2[pt1] + rproducts2[pt1 + 1];
  821.       area2 += tproducts2[pt1 + 2] + tproducts2[pt3 + 2];
  822.       area2 -= rproducts2[pt2 + 1] + rproducts2[pt2];
  823.       area2 -= tproducts2[pt3] + tproducts2[pt1];
  824.       if (area1 != 0.0) area1 = fabs(area2 / area1);
  825.       colors[pt3 + 1] *= area1;
  826.     }
  827.   }
  828.  
  829.   /* Compute the area around the outside of the disk */
  830.   r = rdivisions;
  831.   for (t = 0; t < tdivisions; t++) {
  832.     pt1 = (t ? t-1 : tdivisions-1)*(rdivisions+1) + r - 1;
  833.     pt3 = t*(rdivisions + 1) + r - 1;
  834.     pt2 = ((t+1) % tdivisions)*(rdivisions+1) + r - 1;
  835.     area1 = rproducts1[pt1];
  836.     area1 += tproducts1[pt1 + 1] + tproducts1[pt3 + 1];
  837.     area1 -= rproducts1[pt2];
  838.     area1 -= tproducts1[pt1] + tproducts1[pt3];
  839.     area2 = rproducts2[pt1];
  840.     area2 += tproducts2[pt1 + 1] + tproducts2[pt3 + 1];
  841.     area2 -= rproducts2[pt2];
  842.     area2 -= tproducts2[pt1] + tproducts2[pt3];
  843.     if (area1 != 0.0) area1 = fabs(area2 / area1);
  844.     colors[pt3 + 1] *= area1;
  845.   }
  846.  
  847.  
  848.   delete rproducts1;
  849.   delete tproducts1;
  850.   delete rproducts2;
  851.   delete tproducts2;
  852. }
  853.  
  854. void Unitdisk::fill_trig_tables()
  855. {
  856.   int t;
  857.  
  858.   delete sintable;
  859.   delete costable;
  860.   sintable = new float[tdivisions];
  861.   costable = new float[tdivisions];
  862.   for (t = 0; t < tdivisions; t++) {
  863.     costable[t] = cos(M_2PI * (float)t / (float)tdivisions);
  864.     sintable[t] = sin(M_2PI * (float)t / (float)tdivisions);
  865.   }
  866. }
  867.